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ABSTRACT 

We investigate the impact of massive neutrinos on the distribution of matter in the 
semi-non-Hnear regime (0.1 ^k< 0.6/iMpc^^). We present a suite of large-scale A^- 
body simulations quantifying the scale dependent suppression of the total matter 
power spectrum, resulting from the free-streaming of massive neutrinos out of high- 
density regions. Our simulations show a power suppression of 3.5 — 90 per cent at 
fc~0.6/iMpc~^ for total neutrino mass, Sto^ — 0.05 — 1.9 eV respectively. We also 
discuss the precision levels that future cosmological datasets would have to achieve in 
order to distinguish the normal and inverted neutrino mass hierarchies. 



Subject headings: neutrinos 
scale structure of Universe. 



1 INTRODUCTION 



methods: numerical - large- 



In the standard model of particle physics there are three 
types (flavours) of neutrinos: electron neutrino (z^e), muon 
neutrino (Vf^) and tau neutrino (vt)- Neutrino oscillation ex- 
periments ( [KamLA ND 2008, |SNO|2004[ ) in the past decade 
indicate that at least two neutrino eigentstates have non- 
zero masses. The direct implication of massive neutrinos is 
a non-zero hot dark matter (HDM) contribution to the to- 
tal energy density of the Universe. Being sensitive to the 
mass squared differences between the neutrino eigentstates, 
the oscillation experiments only provide a lower bound 
on the total neutrino mass. Mass splittings of |Am|2| = 
(2.43±0.13) X 10"^ eV^ and Amji = (7.5 9±0.21) x 10"^ eV^ 
( [Adamson et al.|2008||KamLAND|2008[ ) imply a lower limit 
for the sum of the neutrino masses to be 0.05 and 0.1 eV for 



the normal and inverted mass hierarchies (Otten & Wein- 
heimer|["2008 1 , respectively. 



During the radiation era, matter perturbations on the 
sub-horizon scales grow logarithmically. The earlier a mode 
enters the horizon, the more it is suppressed due to the 
decaying gravitational potentials. On the other hand, the 
superhorizon modes do not decay until they enter the hori- 
zon. As a result, the matter power spectrum turns over at a 
scale that corresponds to the one that entered the horizon at 
radiation-matter equality. Neutrinos with mass on the sub- 
eV scale behave as a hot component of the dark matter. Neu- 
trinos stream out of high-density regions into low-density 
regions, thereby damping out small-scale density perturba- 
tions. Massive neutrinos, therefore, suppress the logarithmic 



growth of sub-horizon modes. Extremely low mass neutri- 
nos become non-relativistic after the radiation era is over 
and the free-streaming damping of matter perturbations af- 
fects even those scales that were always outside the horizon 
during the radiation era. 

The redshift-dependent free-streaming comoving wave 
number, kfs, is given by 



Hiz) 



(1) 



where H(z) and iith are the Hubble parameter and the neu- 
trino thermal velocity, respectively. As long as neutrinos 
are relativistic, they travel at the speed of light and their 
free-streaming comoving wave number shrinks at the same 
rate as that of the comoving Hubble wave number (equation 
1). After a neutrino eigentstate becomes non-relativistic, its 
thermal velocity decays as 
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where is the mass of a neutrino eigentstate in e V and the 
present-day photon temperature, T^, is 2.725 K (Komatsu 
eraL][2010l ) 



Thus the free-streaming comoving wave number for 
non-relativistic neutrinos is given by 



; 0.81 



(1 + 2)2 



Vlevy 



/iMpc 



(3) 



For a massive eigentstate, the redshift of non-relativistic 
transition (m^ ~ 3T,^) is given by 



1987 
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(4) 
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After a neutrino eigentstate becomes non-relativistic, kfs be- 
gins to grow as fcfs oc (1 + z)~^^^. Thus, kis passes through 
a minimum, fcnr, which can be shown to be (from equation 
3) 



knr 



0.018 (^)''"(f7.„fe^)^/^Mpc- 



(5) 



For modes with k > kfs, the neutrino density pertur- 
bations are erased. This weakens the gravitational potential 
wells and the growth of cold dark matter (CDM) pertur- 
bations is suppressed. Perturbations are free to grow again 
once their comoving wave numbers fall below fcfs. Modes 
with k < km- are never affected by free-streaming and neu- 
trino perturbations evolve like CDM perturbations. Baryon 
density perturbations, on the other hand, being pressure 
supported, can grow in amplitude only after photon decou- 
pling. At the time of photon decoupling, baryons fall into 
the neutrino-damped dark matter potential wells. Thus, ac- 
curate measurements of the amplitude of clustering of mat- 
ter in the Universe can provide strong upper bounds on the 
mass of neutrinos. 

Section 2 describes how we implement neutrinos in our 
A'^-body simulations. In this section we also discuss the nu- 
merical methods employed in some complementary recent 
studies. In section 3 we discuss the convergence tests for the 
matter power spectrum calculated from our A^-body simula- 
tions. In section 4 we show the impact of massive neutrinos 
on the matter distribution through the total matter power 
spectrum. In section 5 we discuss, based on our A/^-body sim- 
ulations, the precision levels at which future galaxy surveys 
would need to measure the matter power spectrum, in order 
to distinguish between the normal and inverted mass hierar- 
chies. In section 6 we compare our results with the neutrino 
simulations performed by other groups. In section 7 we es- 
timate the errors in our TV-body matter power spectra. We 
present our conclusions in section 8. 



2 IMPLEMENTING NEUTRINOS IN THE 
Af-BODY SIMULATIONS 

Neutrinos in the mass range 0.05 < Em^ < 1 eV have 
present-day free-streaming scales 0.04 < fcfs < 0.3 hMpc~^ 
(150 > Afs > 20/i~^Mpc) and thermal velocities 
3000 > ?;th > 450km/s respectively. Such large ther- 
mal velocities would prevent neutrinos from clustering with 
CDM and baryons, thereby keeping the neutrino perturba- 
tions in the linear regime. As such, in our simulations, we 
have assumed that the non-linear neutrino perturbations 
can be ignored and include the linear neutrino perturbations 
in the initial conditions (ICs) only. 

To generate the ICs for CDM particles and baryons, 
we use the publicly available CA MB code ([Lewis fc Challi- 



2002|| and e nzo1.5 cod^ ( |0'Shea et al.||2004| iNor 



man et al. 20071 - an adaptive mesh refinement (AMR), 



grid-based hybrid code (hydro + TV-Body) designed to sim- 
ulate cosmological structure formation. We use the CAMB 
code to calculate the linear transfer functions for a given 
CDM-|-baryon-|-neutrino-|-A model. The linear density fluc- 
tuation field for CDM particles and baryons is then calcu- 
lated from their transfer functions using ENZOl.5. The initial 
positions and velocities for CDM particles and baryon veloc- 
ities are calculated using the Zel'dovich Approximation (ZA, 
Zel'dovich (19701). Note that we do not have neutrinos in 
our simulations as TV-body particles or as a linear grid. Neu- 
trinos enter our simulations only as neutrino- weighted CDM 
and baryon transfer functions from CAMB. 

The linear matter power spectrum, P^, can be calcu- 
lated at z = as the weighted average of the neutrino (Pt) 
and the combined CDM plus baryon(P^^) linear spectra: 

Pt{k) = [{fc + .MVWk)+f.VPmy : (6) 

where the weights are fi = Sli/fim and = fib + Sic + i^i^- 
The CDM plus baryon power spectrum is 

2 



(7) 



where P^ and P^ are the linear CDM and baryon power 
spectra respectively. The superscript 'L' indicates quantities 
in the linear regime. On smaller scales the matter perturba- 
tions have gone non-linear. So, the non-linear matter power 
spectrum, P^^ , at z — becomes 



(8) 



where. 



PT^ik) = (/c + /b)-' [hVF¥Hk) + .fb/frifc)) ' . (9) 



In equation (8), we calculate P^ at z 

3L 



from TV-body 

simulations and combine it with P^ at 2: = as solved by 
the CAMB code to construct P^{k). Note that we do not 
account for the non-linear neutrino corrections in equation 
(8). Saito et al. ( 2009!) studied the non- linear neutrino per- 



turbations using the higher-order perturbation theory (PT) 
to show that for low neutrino fractions {f^ ^0.05), the am- 
plitude of the non-linear matter power spectrum increases 
by <0.01 per cent at fc~0.2/iMpc"^ at z = 3 and by<0.15 
per cent at fc ~ O.l/iMpc"^ at z = 0. Since at 2 = 0, the 
PT approach to the non-linear matter power spectrum is ex- 
pected to reproduce the TV-body simulation results within 
1 per cent - only for fe^O.l — 0.15/iMpc~'^ (Taruya et al 



20091, the non- linear neutrino corrections at 2 = may 



be somewhat larger on scales we probe in our simulations 
(0.1 < fc < 0.6/iMpc~^) - the estimate of which requires 
multiple particle (CDM-|-baryon-|-neutrino) simulations. 

Numerical studies of the effect of neutrinos on the mat- 
ter distribution have recently been performed independently 
('2OO8); Brandbyge fc Hannestad| p009 



by Brandbyge et al. 



http: / /lea. ucsd.edu/projects/enzo 



2010 1 and Viel et al._ (2010) . Both groups choose similar 
cosmological parameters: (fim = 0.3, fib = 0.05, S7c + S^i^ = 
0.25, = 0.7, ft = 0.7, ria = 1), a 512ft~ ^Mpc box and an 
initial redshift for simulations, Zi = 4 9. [Brandbyge et al.| 
(2008[) and [Brandbyge fc Hannestad] (|2009| |2010[) use a 
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weighted sum of the CDM+baryon transfer functions (since 
they do not have baryons in their simulations) to generate 
ICs for the CDM component using ZA+second-order La- 



grangian perturbation theory (2LPT; Scoccimarro (1998l) 



The Viel et al. (20101 simulations include baryons and use 
ZA to generate ICs. Both groups include neutrinos in their 
iV-body simulations either as A''-body particles, as a linear 
grid or use a hybrid method where neutrinos are treated 
as grid or particles depending on their thermal motion. In 
the grid-based implementation, the neutrino grid is evolved 
linearly and does not include the non-linear corrections. 
The particle-based implementation accounts for the non- 
linearities by including the coupling between the gravita- 
tional potential and neutrinos. 

Brandbyge & Hannestad (20091 (their fig. 1, middle 



panel) show that the error from neglecting non-linear neu- 
trino perturbations at z = is at most 1.25 per cent level 
at A; ~ 0.25/iMpc"^ for T,m„ = 0.6 eV. Also, the error be- 
tween the grid and particle representations is shown to be- 
come smaller on small scales. Specifically, the two repre- 
sentations converge for fc ^ 0.2 /iMpc~^. This is attributed 
to the fact that the neutrino white noise (due to the fi- 
nite number of neutrino A''-body particles) contribution to 
the matter power spectrum dominates only on ever smaller 
scales as the CDM perturbations grow at low redshifts. |Viel| 
et al. (20101 (their fig. 2, right panel) show that the non- 
linear correction at z = may be as high as 6 per cent 
at ~ l/iMpc~^ for Em^ = 0.6 eV and the agreement be- 
tween the grid and particle representations begins to im- 
prove only at fc ^ 1 /iMpc~^. The discrepancies between the 
results from the two groups worsens significantly when the 
above comparison is done at z = 1. These large discrepancies 
can not be solely due to the absence/presence of baryons or 
whether ZA or ZA-I-2LPT is used to generate the ICs since 
(i) the baryons closely trace the CDM distribution on scales 
A: ^ 1 /iMpc"^ and (ii) ZA or ZA+2LPT do not affect the fi- 
nal results significantly when the simulations start at a high 
redshift {zi = 49). The extent and the scale-dependence of 
non-linear neutrino corrections are still being researched. 



3 iV-BODY SIMULATIONS: OPTIMIZING 
BOXSIZE AND NUMBER OF PARTICLES 

We performed TV-body simulations with the ENZOl.5 code. 
The code allows us to choose the geometry (box size, num- 
ber of particles), the cosmology (ilm, f2b, SIa, fii/), the am- 
plitude of fluctuation on 8 Mpc scale: ag, the primordial 
spectral index: ris and the initial redshift: Zi. We kept AMR 
off (no adaptive mesh refinement) since it does not signif- 
icantly affect the scales of interest. Throughout this paper 
we assume the 7-yr Wilkinson Microwave Anisotropy Probe 



[WMAP; Larson et al. (2010l) central parameters: ilm = 
0.266, fib = 0.044, = 0.734, h = 0.71 and = 0.963 for 
the matter, baryonic and cosmological constant normalized 
densities, the Hubble constant and the primordial spectral 
index respectively. We vary Jl^ such that Qcdm + ^v ~ 0.222. 
The simulation parameters are listed in Table 1. In order 
to suppress sampling variance of the estimated power spec- 
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Table 1. Simulation parameters. All simulations were started at 
a redshift of Zi = 20 and stopped at z = 0. We ran eight indepen- 
dent simulations for each row to suppress sampling variance. 



trum, for each row we ran eight simulations by changing the 
seed to generate the ICs. 



First, we had to select an appropriate geometry (box 
size and the number of CDM/gas particles) for which the 
matter power spectrum converges to 1 per cent accuracy 
in the semi-non-linear regime (0.1 ^ k <, 0.6/iMpc~^). 
The largest mode that can fit in a 200/i~^Mpc box IS 
k ~ 0.03/iMpc~^ and the matter power spectrum is suffi- 
ciently linear on these scales. One can choose bigger volumes 
but unless the number of particles is also increased accord- 
ingly, it leads to a poor mass resolution. Also, A^-body simu- 
lations suffer from a discreteness problem that arises due to 
the finite number of macroparticles used to sample the mat- 
ter distribution in the universe. Thus, given any theoretical 
cosmological model, the ICs are always undersampled. 

The smallest scale for which the power spectrum can be 
resolved accurately is related to the Nyquist wavenumber, 
fcNy, given by: 



AN, 



part J 



xl/3 



(10) 



Given a combination of the number of particles and the box 
size, the power spectrum is dominated by shot noise for k ^ 
fcNy. For A'^cdm = 64^ particles in a 200/i~^Mpc box, fc^y 
is l.Ol/iMpc"^, while the semi-non-linear modes of interest 
are 0.1 ^ fc ^ 0.6/iMpc"\ Thus iVcdm = 64^ particles in a 
200ft~^Mpc box seems a reasonable combination to start 
with. 

The number of gas particles fixes the root grid that de- 
termines the force resolution for the simulation, enzo uses a 
particle mesh technique to calculate the gravitational po- 
tential on the root grid (O'Shea et al. 2005). Forces are 



first computed on the mesh by finite-differencing the gravi- 
tational potential and then interpolated to the dark matter 
particle positions to update the particle's position and ve- 
locity information. This methodology requires that the root 
grid be at least twice as fine as the mean interparticle separa- 
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k (h Mpc-') 

Figure 1. Matter power spectrum at z = for undersam- 
pled ICs at Zi = 20 with Afcdm = 64^ - solid (red), 128^ - 
long dash-dotted (green), 256^ — dashed (blue) and 512"^ — 
long-dashed (cyan). The vertical lines are the fe^y wavenum- 
bers for 64^, 128^, 256^ and 512^ CDM particles. Also plotted 
(dot— dashed line) is the linear theoretical power spectrum. For 
k > /cNy, particle shot noise dominates the true power spectrum. 




0.1 1 

k (h Mpc-i) 

Figure 2. Same as Fig. [l] expressed as fractional suppression of 
the matter power spectrum at z = when 64'^ — solid (red) , 128^ — 
long dash-dotted (green) and 256^ — dashed (blue) CDM par- 
ticles are used to sample the ICs w.r.t the case where 512^ — 
long-dashed (cyan) CDM particles are used. Qi, = for all four 
cases. The error bars correspond to eight simulations with differ- 
ent seeds for the ICs. 



tion to obtain accurate forces down to the scale of the mean 
interparticle spacing. A coarse root grid renders the forces 
on the scale of the mean interparticle spacing, inaccurate. 
This explains our choice of A'^gas = 512'^. 

Fig. [Tlshows the matter power spectrum at z = when 
Afcdm = 64^, 128^, 256^ and 512^ particles are used to sam- 
ple the ICs (f2i/ = for all four cases). Beyond the Nyquist 
wavenumbers, represented by vertical lines [64"^ - solid (red), 
128'^ - long dash-dotted (green), 256'^ - dashed (blue) and 
512'^ - long-dashed (cyan)], the power spectra become in- 
creasingly inaccurate due to particle shot noise contribution. 
Fig. [2] shows the fractional suppression of the matter power 
spectrum at z — 0. For k <, 0.6/iMpc~^, the error due to 
undersampling the ICs is ^5 per cent for the 64^ run, ^0.5 
per cent for the 128'' run and negligibly small for the 256^ 
run. To keep the undersampling error at fe = 0.6 /iMpc^^ 
below 0.5 per cent, we narrowed down to a combination of 
A^cdm = 256^, A^gaa = 512^ in a 200/i"^Mpc box to investigate 
the effect of massive neutrinos on the matter power spec- 
trum in the semi-non-linear regime (0.1 < k < 0.6^Mpc~^). 
Finally, we checked the smallest scales that are accurately 
resolved by the 200fe~^Mpc box. Towards this, we ran 
eight simulations in a 100/i~^Mpc box with A'^cdm = 256"', 
Afgas = 512^. In Fig. [3] we plot the power spectrum from 100 
and 200^~^Mpc boxes. The matter power spectrum from 
100/i~^Mpc box simulations begins to show excess power 
for fc^l/iMpc"^ The non-linear evolution of perturbations 
on scales k ^ 1/iMpc^^ is missed in the 200/i~^Mpc box 
simulations. The spectrum from 200 h~^Mpc box simula- 
tions show convergence at 1 per cent level for k<,l hMpc~^ 
(Fig.|4f. 



4 IMPACT OF MASSIVE NEUTRINOS 

The contribution of massive neutrinos to the present-day 
critical energy density is given by: 



(11) 



94.22/i2 

where Em^ is the sum of the masses of all neutrino 
eigentstates. In this section we consider four neutrino 
models: Q,^ — 0, 0.01, 0.02 and 0.04 corresponding to 
Errii^ = 0, 0.475, 0.95 and 1.9 eV, respectively. We as- 
sume three degenerate neutrino eigentstates, so that 

In Fig. [5] we show slices of the baryon density field at 
2 = extracted from 200/i"^Mpc box with A''cdm = 256^, 
A'^gas ~ 512^. The top panel is from a simulation without 
neutrinos, the middle and the bottom panels correspond to 
simulations with $1,^ = 0.02 and 0.04 respectively. All slices 
are 200 Mpc wide. The slices show the baryonic mass 
averaged over the volume of a grid cell. Each grid cell in our 
simulations is ~391 h~^kpc. 

As neutrinos become more massive, the suppression 
in the growth of density perturbations becomes clear by 
the relatively diffused density filaments. The baryon den- 
sity fields in the middle and the bottom panels are less 
evolved relative to the massless neutrino (top panel) case. 
The gravitational potential wells are much deeper in the 
top panel. This is evident from the voids (dark blue regions) 
which are more underdense in the top panel compared to 
the voids in the lower panels. To quantify the difference 
between simulations with and without massive neutrinos. 
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Figure 3. Matter power spectrum at 2 = from 
100 fc^^Mpc - solid (green) and 200 /i-^Mpc - dashed (blue) 
box simulations. The linear theory spectrum (dot— dashed line) is 
also shown. The vertical dashed line is the maximum wavenumber 
up to which the power spectrum from 200/i~^Mpc box simula- 
tions can be trusted at 1 per cent level. 




k (h Mpc-i) 

Figure 4. Same as Fig. |3] expressed as fractional suppres- 
sion of the matter power spectrum at 2: = as a function of 
the box size. Spectrum from 100/i~^Mpc — solid (green) and 
200/i~^Mpc — dashed (blue) box agree at 1 per cent level for 
A:^l/iMpc-i. 



we measure the total matter power spectrum by converting 
the positions of the CDM and gas particles into 512'^-point 
grids of densities using a Cloud-In-Cell (CIC) interpolation 
scheme. We do not compensate for the smoothing effect in- 
troduced by the CIC filtering since the smoothing affects 
scales that are close to the Nyquist wavenumber which for 



our choice of parameters (A'^gas = 512^, Box=200 ft~^Mpc) 
is ^Ny ~ 8.04/iMpc~^, while the semi-non-linear modes of 
interest are 0.1 ^ fc^O.6 /iMpc~^. The density fields are fast 
Fourier transformed to calculate P^^{k) and P^^{k) - the 
non-linear power spectrum for baryons and CDM respec- 
tively. We then construct P^^{k) at z = using equations 
(8) and (9). To suppress sampling variance of the estimated 
P{k), we take the average P{k) from eight simulations. 

Fig. [6] shows the matter power spectrum at z = from 
simulations and linear theory (dot -dashed lines) as a func- 
tion of neutrino mass for the four neutrino models: Sl^ = 
0(Em^ = OeV) - solid (red), = 0.01 (Em„ = 0.475 eV) 
- long dash-dotted (green), fl^ = 0.02 (Em„ = 0.95 eV) - 
dashed (blue) and fi;^ = 0.04 (Em^ = 1.9 eV) - long-dashed 
(cyan). The simulation spectra are significantly above the 
linear theory predictions at high k. The linear theory pre- 
dictions break down for k ^ 0.1/iMpc~^ (A ^ 60 /I'^Mpc). 
Also, as the total neutrino mass is increased (keeping the 
number of degenerate neutrino eigentstates fixed at three), 
the matter power spectrum is further suppressed. Since neu- 
trino eigentstates with higher mass constitute a larger frac- 
tion of the total energy density, they are more effective in 
damping small-scale power than low mass neutrinos. 

In Fig. [7] we plot the fractional difference between the 
matter power spectra with and without massive neutrinos, 
from the simulations as well as the linear theory predictions. 
The linetypes for the spectra are the same as in Fig. [6] The 
linear theory predicts a nearly scale-independent suppres- 
sion for A: ^ 0.2 ftMpc"^. On the other other hand, the non- 
linear power spectra from the simulations show an enhanced 
suppression for fe 0.1 /iMpc~^. At fc ~ 1 /iMpc"'^, the non- 
linear spectra are ~ 10 per cent more suppressed compared 
to the corresponding linear spectra. 



RESOLVING NEUTRINO MASS 
HIERARCHY FROM SIMULATIONS 

The mass splittings of lAmial = (2.43 ± 0.13) x 10"^ eV 



and Amgi = (7.59 ± 0.21) x lO"'^ eV^ (Adamson et al. 



2008 KamLAND 20081 allow for two possible neutrino 
mass hierarchies: normal {m^ > m2 > mi) and inverted 
(m2 > mi > ma). For Em„ > 0.4 — 0.5 eV, all neu- 
trino eigentstates are essentialy degenerate, the mass of each 
eigentstate being rrii, « Em;^/3. However, for smaller Em^, 
the individual eigentstate masses differ significantly in the 
normal and inverted hierarchies. The free-streaming comov- 
ing wave number, fcnr, is a function of the mass of each neu- 
trino eigentstate (see equations 4 and 5). As the mass is 
increased, it becomes non-relativistic earlier and the free- 
streaming scale gets shorter. The mass dependence of k^ 
means that the matter power spectrum is modified differ- 
ently for eigentstates with different masses. This makes the 
matter power spectrum a powerful tool to distinguish be- 
tween the normal and inverted hierarchies. In this section 
we discuss the precision levels above which the power spec- 
trum from future galaxy surveys should be able to resolve 
between the two mass hierarchies. 

The mass splittings of jAmial = (2.43±0.13) x 10"^ eV^ 
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I 



I 



Figure 5. Slices of baryon density distribution. All slices are 
200/i~^Mpc wide and show the baryonic mass averaged over the 
volume of a grid cell. Each grid cell is ~ 391 /i~^kpc. The top panel 
shows a simulation without neutrinos. The middle and the bot- 
tom panels are taken from simulations with f2^ = 0.02 (Sm^ = 
0.95eV) and 0„ = 0.04 (Sm^ = 1.9eV). The baryon density 
fields in the middle and the bottom panels are less evolved rel- 
ative to the no-neutrino (top panel) case. The simulations were 
run with A^cdm = 256^, Afgas = 512'^. The density projections were 
made using YT: an analysis and visualization tool l |Turk|2008[ l . 




0.1 1 

k (h Mpc-i) 

Figure 6. Matter power spectrum at z = from simulations 
and linear theory (dot-dashed lines) as a function of neutrino 
mass. The four neutrino models are: fi^ = 0(Emiy = OeV) - 
solid (red), Q.^ = 0.01 (Sm„ = 0.475 eV) - long dash-dotted 
(green), Q.^ = 0.02 (Sm^ = 0.95 eV) - dashed (blue) and 0.u = 
0.04 (Srrti, = 1.9 eV) - long-dashed (cyan). The vertical dashed 
line is the maximum wavenumber up to which the power spectra 
from 200/i~^Mpc box simulations are valid at 1 per cent level. 
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0.1 
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Figure 7. Fractional difference between the matter power spec- 
tra with and without massive neutrinos at z = 0, from the simu- 
lations and the linear theory predictions (dot-dashed lines). The 
four neutrino models are: = 0(Sm^ = OeV) - solid (red), 
= 0.01 (Sm^ = 0.475 eV) - long dash-dotted (green), Q.^ = 
0.02 (Sm,. = 0.95 eV) - dashed (blue) and = 0.04 (Sm^ = 
1.9 eV) - long-dashed (cyan). The error bars correspond to eight 
simulations with different seeds for the ICs. 
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and Amil = (7.59 ± 0.21) x 10"^ eV^ imply that the lower 
bounds on the total neutrino mass are Em^ — 0.05 and 
0.1 eV for the normal and inverted mass hierarchies respec- 
tively. We performed A'^-body simulations for Em^ = 0.05 
and 0.1 eV. For Sm^ = 0.05 eV, we assumed 1 massive and 
2 massless eigentstates (mimicking the normal hierarchy). 
For Em^ — 0.1 eV, we assumed 2 massive and 1 massless 
eigentstate (mimicking the inverted hierarchy). In Fig. [s] 
we show the fractional suppression in the power spectrum 
for two neutrino models: Q^, — 0.001 (Emi/ = 0.05 eV) - 
long dash-dotted (green) and 0,^ = 0.002 (Emj, — 0.1 eV) 
- dashed (blue). The growth of structure formation is sup- 
pressed by as much as 3.5 per cent (7.5 per cent) at ~ 
0.6/iMpc~^ for the two models. The measurement errors in 
the power spectrum from future galaxy surveys are expected 
to be at the 1 per cent level. In case future surveys constrain 
Em^ < 0.1 eV with sufficient precision, that would rule out 
the inverted mass hierarchy. The current constraint from the 



7-yr WMAP data alone (Larson et al. 20101 is Em^ < 1.3 eV 



(95 per cent CL). At this level, it is not possible to discrim- 
inate between the normal and inverted hierarchies since all 
eigentstates are essentially degenerate. 

Next, we consider a scenario with Em^ — 0.1 eV, at 
which the difference between the normal and inverted hier- 
archies is most prominent. We ran A'^-body simulations in the 
following three ways: (i) (A^massive = 3, A'^dogcn = 3) where 
A'massive is the uumber of massive eigentstates and A^degen 
is the degeneracy amongst the massive eigentstates. This 
combination corresponds to rrii, = TjUIv/?, = 0.033 eV; (ii) 
(A^niassivc = 2, A'^dcgon = 2), this is the inverted hierarchy sce- 
nario with one massless and two equally massive eigentstates 
(m, ~ 0.05, 0.05, OeV); (iii) (Ar^naBBivc = 3, A^dcgcn — 2), 
this is the normal hierarchy scenario with three massive 
eigentstates (m^ ~ 0.056, 0.022, 0.022 eV). Note that case 
(i) is meaningless at Em,^ = 0.1 eV given that |Am32| = 
(2.43±0.13)xlO"^eV^ and Amii = (7.59±0.21) x 10"'^ eV^ 
We include case (i) for illustrative purposes only. 

In Fig. [9] we plot the matter power spectrum for 
cases (i), (ii) and (iii) divided by the spectrum for case 
(i). The linear theory predictions are shown by dot-dashed 
lines. Since non-linearities become important only for k ^ 
O.l/iMpc"^, we have plotted the theoretical power spec- 
trum for k < 0.1 /iMpc~^, calculated using the CAMB code. 
The suppression from simulations is ~ 0.05 — 0.2 per cent 
higher than the linear predictions. The inverted hierarchy 
- dashed line (green) shows excess power for wavenumbers 
0.001 < k < 0.02 ftMpc"^ and an enhanced suppression of 
~0.5 per cent at fc~l/iMpc~^ relative to case (i). This can 
be explained by the fact that in case (ii) Emi, = 0.1 eV is 
shared equally between two eigentstates, while in case (i) 
Em^ = 0.1 eV is shared equally between three eigentstates. 
Each eigentstate is more massive in case (ii), thereby making 
the free-streaming length shorter compared to that in case 
(i). Higher mass neutrinos are better at wiping out small- 
scale perturbations and their shorter free-streaming length 
implies that the spatial extent of damping is limited. 

Another factor contributing to the appearance of Fig. [9] 
is a shift in the radiation-matter equality redshift. Higher 
mass neutrinos become non-relativistic at higher redshifts 
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Figure 8. Same as Fig. [7] but for neutrino models with much 
lower neutrino mass: Cj, = 0.001 (Sm^ = 0.05 eV) - long dash- 
dotted (green) and = 0.002 (Sm^ = 0.1 eV) - dashed (blue). 



and start contributing to before low mass neutrinos do. 
This shifts the radiation-matter equality epoch to a higher 
redshift and reduces the scale corresponding to the one that 
entered the horizon at radiation-matter equality. The modes 
entering the horizon after radiation-matter equality grow 
linearly (as opposed to logarithmically during the radia- 
tion era) which contributes to the excess power [compare 
dashed (green) and solid (red) lines in Fig. [9] for wavenum- 
bers 0.001 <k < 0.02/iMpc~^. The same reasoning can be 
applied to the normal hierarchy - long dash-dotted line 
(blue). At Em^ — 0.1 eV, precision better than 0.5 per cent 
would be needed in measuring the matter power spectrum 
to discriminate between the normal and inverted hierarchies. 
For Em„ > 0.2 eV all eigentstates become degenerate, this 
would make it extremely difficult for a future survey to re- 
solve the two hierarchies. 



6 COMPARISON WITH RECENT 
NUMERICAL STUDIES 

In this section we compare the estimated overall suppres- 
sion of the matter power spectrum due to massive neutrinos 
from our Af-body simulations with the results obtained by 
[Brandbyge et aLlpOOS) ) and |Viel et"aL] ( |2010[ ). In linear the- 
ory, the suppression of the matter power spectrum ampli- 
tude is approximately given by AP/P ~ —8/^ (Hu, Eisen- 



stein, & Tegmark 19981. Numerical simulations, however, 

inced in the non- 
we plot the frac- 



show that the neutrino suppression is enhanced in the non- 
linear regime (fc ^ 0.1 feMpc"^) 



10 



In Fig. 

tional difference between the matter power spectra with and 
without massive neutrinos at 2: = 0, from the simulations as 
well as the linear theory predictions (dash-dotted lines) for 
four neutrino models: Sl^ — 0.001 (Errii/ — 0.05 eV) - dot- 
ted (green), Sl^ = 0.002 (Em^ = 0.1 eV) - dashed (blue). 
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Figure 9. Matter power spectrum for normal - long dash-dotted 
line (blue) and inverted - dashed line (green) hierarchies divided 
by the matter power spectrum for mj, = - solid line 

(red). The linear theory predictions are shown by dot-dashed 
lines. The neutrino model considered here is Em^ = OeV. The 
individual masses for the three eigentstates are (m„ ~ 0.05, 0.05 
and OeV) for the inverted hierarchy and (m,^ ^ 0.056,0.022 and 
0.022 eV) for the normal hierarchy. The inverted hierarchy shows 
more damping of small-scale power than the normal hierarchy. 



Figure 10. Fractional difference between the matter power spec- 
tra with and without massive neutrinos at z = 0, from the sim- 
ulations and the linear theory predictions (dash-dotted lines). 
The four neutrino models are: = 0.001 (Sm^ = 0.05 eV) - 
dotted (green), = 0.002 (Emi, = 0.1 eV) - dashed (blue), 
= 0.01 (Sm^ = 0.475 eV) - long-dashed (cyan) and = 
0.02 (Smiy = 0.95 eV) — long dash-dotted (magenta). The maxi- 
mum relative suppression of AP/P~ — lO/i^ is shown as short hor- 
izontal dotted lines. The horizontal (red) dotted line for Em,/ = 
0.95 eV is at AP/P~-8.6/^. 



n^^O.Ol (Em,!/ = 0.475 eV) - long-dashed (cyan) and n„ = 
0.02 (Em^, = 0.95 eV) - long dash-dotted (magenta). We 
found a maximum non-linear suppression of AP/P~ — lO/i/ 
for neutrino masses Em^ = 0.05, 0.1, 0.475 eV. Although 
we ran our simulations with a slightly different set of cos- 
mological parameters, Brandbyge et al. ( 2008 1 measured 
AP/P~-9.8/^ for Em^ < 0.6 eV while |Viel etaT] poTo I 
reported AP/P ~ -9.5/^ at z = 0. For Em,. = 0.95 eV, 
we get AP/P ~ -8.6/1, while |Viel et al.| ( [Mof reported 
AP/P~ -8/^ for Emj. = 1.2 eV. The scale at which the sup- 
pression turns over, fcnr, moves from fcnr~0.6 — 0.7/iMpc~^ 
for Em,^ = 0.05 eV to fc^r ~ l/iMpc^^ for Em^ = 0.95 eV. 
The turnover may be related to the non-linear collapse of 
structures as discuss^ 
ported fcnr ~ 1 /iMpc" 



structures as discussed in Brandbyge et al. ( 2008 1 who re 
_i 



7 MATTER POWER SPECTRUM ERROR 
ESTIMATES 

In our A^-body simulations, we have implemented neutrinos 
in the ICs only. Neutrino-weighted CDM and baryon trans- 
fer functions from camb were used to generate the ICs for 
CDM particles and baryons. To construct P^iji) at 2: = 0, 
we used equation (8), where Pf^ at z = from A-body 
simulations was combined with P^ at z = as solved by 
the CAMB code. This methodology introduces errors in the 
estimated matter power spectrum for two reasons: (i) the 
linear neutrino perturbations were taken into account only 
at the initial (zi = 20) and the final (2 = 0) redshifts. 



There is no feedback from the neutrinos on to the CDM 
component in our A-body simulations, (ii) the non-linear 
evolution of neutrino perturbations was not accounted for 
in our A-body simulations. While the extent of non-linear 
neutrino corrections to the matter power spectrum is still be- 



ing studied, we use Brandbyge et al. ( 2008 1 and Brandbyge 



& Hannestad ( 2009 1 to estimate the errors in our A-body 
spectra. [Brandbyge fc Hanncstad ( 2009 ) describe the linear 
neutrino density on a grid and evolve this density forward in 
time using linear theory. The neutrino contribution is added 
to the CDM component when calculating the gravitational 
forces. Thus, the linear neutrino component is accounted 
for recursively over the redshift range over which the matter 
power spectrum is to be evolved. [Brandbyge et al.] ( |2008| ) 
(their fig. 7, left panel) show that the matter power spec- 
trum is underevolved by ~ 3 per cent for Emj, < 0.6 eV 
on scales k > 0.2 /iMpc~^ when the neutrino grid is ne- 
glected. Accordingly, our matter power spectrum estimates 
are expected to be underevolved by roughly Si 5, 3 and 0.1 
per cent for Em^ = 0.95, 0.475 and 0.1 eV, respectively, for 



k ^0.2/iMpc at 2 = 0. Fig. 1 in Brandbyge & Hannes 



tad ( 2009 1 shows that the power is further suppressed by 
~ 5 per cent for Em^ < 1.2 eV at fc 0.2 - 0.3/iMpc"^ 
when the neutrino non-linearities are neglected. Overall, we 
estimate our A-body spectrum errors to be 10,4 and 0.1 
per cent for Em^ = 0.95,0.475 and 0.1 eV, respectively, for 
fe^0.2/iMpc"^ at 2 = 0. 
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8 DISCUSSION AND CONCLUSIONS 



9 ACKNOWLEDGMENTS 



In this paper we simulated the matter power spectrum at 
« = in order to study how massive neutrinos impact struc- 
ture formation. The most important factors in obtaining an 
accurate power spectrum are (i) the Nyquist wavenumber, 
which depends on the simulation box size and the number 
of particles and (ii) the force resolution, which depends on 
the size of the root grid. Above the Nyquist wavenumber, 
the power spectrum is dominated by shot noise. For the 
semi-non-linear modes (0.1 ^ k <, 0.6/iMpc~^), we found 
that A'cdm = 256^ in a 200/i~^Mpc box is enough to keep 
the sampling errors under 0.5 per cent. We used a root grid 



of N^. 



512 , which is twice as fine as A'cdm, to accu- 



rately calculate the gravitational forces down to the scale 
of the mean interparticle spacing. We also found that the 
non-linear evolution of perturbations are accurate to within 
1 per cent level only for the scales ^ 1 /iMpc"^ when using 
200/i~^Mpc box. Probing smaller scales with higher preci- 
sion requires a smaller simulation box or a finer root grid. 

We have presented a suite of iV-body simulations show- 
ing the effect of massive neutrinos in the range Sl^ = 
0.001 - 0.04 (Em^ = 0.05 - 1.9 eV) on the distribution 
of matter. Massive neutrinos smooth the neutrino density 
field on sub-free-streaming scales. This makes the gravita- 
tional potential wells shallower than their counterparts in 
a pure ACDM universe, leading to a suppressed growth of 
structure formation. The power is suppressed by as much as 
3.5-90 per cent at fc~0.6/iMpc"^ for Sm,^ = 0.05- 1.9eV 
respectively. In our simulations, we include neutrinos as 
neutrino-weighted CDM and baryon transfer functions at 
the starting redshift, Zi = 20. We have neglected the non- 
linear neutrino corrections to the matter power spectrum 
which may be as high as 1.25 per cent for Em^ = 0.6 eV 
and 5 per cent for Emi/ = 1.2 eV as measured by |Brand^ 



byge & Hannestad ( 2009 1 . Although direct comparison of 



our A'^-body results with those from [Brandbyge fc Hannes-| 
tadl pOOgI was not possible since we ran our simulations 



with a slightly different set of cosmological parameters and 
TiTTiv, nevertheless, we expect our A^-body power spectra to 
be in error by ^ 10, 4 and 0.1 per cent for Em^, = 0.95, 0.475 
and 0.1 eV, respectively, for fc ^ 0.2 /iMpc~^ a,t z — Q. We 
found an overall suppression of power from our simulations 
at z = to be AP/P~ -10/^ for Sm^ < 0.5 eV which is 



slightly higher than the results of Brandbyge et al. ( 2008 1 



and |Viel et al.| ( [20T0| ) who reported AP/P ~ -9.8/^ and 
AP/P~-9.5/^ respectively for Em^ < 0.6 eV. 

As part of the Sloan Digital Sky Survey-Ill, the Baryon 



Oscillation Spectroscopic Survey (BOSS; Ross et al. 2010 1 
is expected to measure the power spectrum with precisions 
at which fi^ ~ 0.01 {T.m^ ~ 0.475 eV) could be ruled out. 
This would significantly improve the current 7-yr WMAP 
data alone constraint of Em^ < 1.3 eV. If Em^ constraints 
from cosmology get as low as 0.1 — 0.2 eV, it will open up a 
possibility to resolve the normal and inverted mass hierar- 
chies, though the matter power spectrum would need to be 
determined with precision levels well below 0.5 per cent. 



Computations described in this work were performed us- 
ing the ENZO code developed by the Laboratory for Com- 
putational Astrophysics at the University of California in 
San Diego ( http: / /lca.ucsd.edu|. We thank the users of yt 
(python-based package for analysing enzo datasets) and 
ENZO for useful discussions and guidance towards running 
and analyzing simulations. We thank the referee for a useful 
and constructive report. This work was supported by the 
National Science Foundation through TeraGrid resources 
provided by the NCSA and by a grant from the Research 
Corporation. HAF has been supported in part by an NSF 
grant AST-0807326, by the University of Kansas General 
Research Fund (KUGRF) and acknowledges the hospitality 
of University College, London and Imperial College in the 
UK and the Institut dAstrophysique de Paris, France. 



REFERENCES 

Adamson P., Andreopoulos C, Arms K. E., Armstrong R., 
Auty D. J., Ayres D. S., Bailer B., Barnes P. D., Barr G., 
Barrett W. L., Becker B. R., Bellas A., Bernstein R. H., 
Bhattacharya D., Bishai M., Blake A., Bock G. J., Boehm 
J., Boehnlein D. J., Bogert D., Bower C, Buckley-Geer 
E., Cavanaugh S., Chapman J. D., Cherdack D., Childress 
S., Choudhary B. C, 2008, Phys. Rev. Lett., 101, 131802 

Brandbyge J., Hannestad S., 2009, Journal of Cosmology 
and Astro-Particle Physics, 5, 2 

— , 2010, Journal of Cosmology and Astro-Particle Physics, 
1, 21 

Brandbyge J., Hannestad S., Haugb0lle T., Thomsen B., 
2008, Journal of Cosmology and Astro-Particle Physics, 
8, 20 

Hu W., Eisenstein D. J., Tegmark M., 1998, Physical Re- 
view Letters, 80, 5255 

KamLAND, 2008, Physical Review Letters, 100, 221803 

Komatsu E., Smith K. M., Dunkley J., Bennett C. L., Gold 
B., Hinshaw G., Jarosik N., Larson D., Nolta M. R., Page 
L., Spergel D. N., Halpern M., HiU R. S., Kogut A., Limon 
M., Meyer S. S., Odegard N., Tucker G. S., Weiland J. L., 
Wollack E., Wright E. L., 2010, ArXiv e-prints 

Larson D., Dunkley J., Hinshaw G., Komatsu E., Nolta 
M. R., Bennett C. L., Gold B., Halpern M., HiU R. S., 
Jarosik N., Kogut A., Limon M., Meyer S. S., Odegard 
N., Page L., Smith K. M., Spergel D. N., Tucker G. S., 
Weiland J. L., WoUack E., Wright E. L., 2010, ArXiv e- 
prints 

Lewis A., Chalhnor A., 2002, Phys. Rev. D, 66, 023531 
Norman M. L., Bryan G. L., Harkness R., Bordner J., 

Reynolds D., O'Shea B., Wagner R., 2007, ArXiv e-prints 
O'Shea B. W., Bryan G., Bordner J., Norman M. L., Abel 

T., Harkness R., Kritsuk A., 2004, ArXiv Astrophysics 

e-prints 

O'Shea B. W., Nagamine K., Springel V., Hernquist L., 

Norman M. L., 2005, ApJS, 160, 1 
Often E. W., Weinheimer C, 2008, Reports on Progress in 

Physics, 71, 086201 
Ross N., Sheldon E. S., Myers A. D., Yeche C, Richards 



© 0000 RAS, MNRAS 000, 000-000 



10 Agarwal & Feldman 



G. T., McMahon R. G., Hennawi J. F., Lee K., Wood- 
Vasey W. M., Woyant A., Pctitjcan P., Eisenstein D. J., 
Nichol R. C, Padmanabhan N., Schlegel D. J., Schneider 
D. P., Strauss M. A., Weinberg D. H., White M., 2010, 
in Bulletin of the American Astronomical Society, Vol. 41, 
Bulletin of the American Astronomical Society, pp. 517 — h 

Saito S., Takada M., Taruya A., 2009, PRD, 80, 083528 

Scoccimarro R., 1998, MNRAS, 299, 1097 

SNO, 2004, Physical Review Letters, 92, 181301 

Taruya A., Nishimichi T., Saito S., Hiramatsu T., 2009, 
PRD, 80, 123503 

Turk M., 2008, in Proceedings of the 7th Python in Science 
Conference, Varoquaux G., Vaught T., Millman J., eds., 
Pasadena, CA USA, pp. 46 - 50 

Viel M., Hachnelt M. G., Springcl V., 2010, Journal of Cos- 
mology and Astro-Particle Physics, 6, 15 

Zel'dovich Y. B., 1970, aap, 5, 84 



© 0000 RAS, MNRAS 000, 000-000 



